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ABSTRACT 


A method  for  computing  three-dimensional  flow  over  an  ogival  body  at  angle 
of  attack  is  described.  An  approximate  set  of  governing  equations  is  given  for 
viscous  flows  which  have  a primary  flow  direction.  A two-level  second-order 
accurate  marching  procedure  is  presented  for  general  equations.  With  this  procedure, 
a three-dimensional  turbulent  flow  can  be  solved  in  any  coordinate  system  by  marching 
along  the  assumed  primary  flow  direction.  General  tube-like  coordinates  are 
developed  for  a class  of  geometries  applicable  to  flows  between  tubular  surfaces. 

The  coordinates  are  then  specialized  to  the  flow  field  bounded  between  an  ogival  body 
at  angle  of  attack  and  its  bow  shock.  Unlike  the  ogival  body  surface,  the  bow  shock 
surface  is  not  known  in  advance  of  the  solution  but  instead  must  be  computed  as  the 
solution  develops.  One  marching  step  of  the  solution  process  is  broken  down  into 
several  steps.  First,  the  bow  shock  surface  is  discretely  extended  by  an  iteration 
of  explicit  local  solutions.  The  bow  shock  surface  is  then  smoothly  extended  to 
provide  a best  fit  to  the  discrete  shock  data.  Tube-like  coordinates  are  generated 
and  finally  the  second  order  numerical  scheme  is  applied  to  advance  the  fully  viscous 
solution  to  the  next  station.  In  addition,  some  preliminary  computational  results 
were  obtained.  Specifically,  the  code  was  applied  to  subsonic  boundary  layers, 
purely  supersonic  flow  with  shocks,  and  mixed  subsonic-supersonic  flow. 
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A Method  for  Computing  Flows  Over  an  Ogival  Body 

INTRODUCTION 


The  accurate  prediction  of  both  the  pressure  distribution  and  heat  transfer 
loads  about  an  ogival  body  is  an  important  consideration  in  the  design  of  super- 
sonic missiles.  The  combination  of  inviscid  flow  theory  and  three-dimensional 
boundary  layer  theory  may  be  adequate  to  predict  many  flows  about  ogival  bodies  at 
small  angles  of  attack;  however,  the  separate  use  of  these  analyses  are  inadequate 
at  larger  angles.  In  such  cases,  a strong  viscous-inviscid  interaction  occurs  on 
the  lee  side  of  the  body;  this  leads  to  the  formation  of  a pair  of  vortices 
symmetric  about  the  lee  body  generator.  But  then,  under  these  conditions,  an 
accurate  flow  field  prediction  requires  an  analysis  of  the  mutual  interaction  between 
the  viscous  layer  and  the  nominally  inviscid  flow.  The  formation  of  vortices  will 
often  lead  to  an  increase  in  the  body  lift  to  drag  ratio  and  in  local  heat  transfer 
rate.  An  accurate  method  of  predicting  the  interacting  flow  field  here  is  necessary 
to  insure  both  the  effective  operation  and  structural  integrity  of  supersonic 
vehicles. 

Successful  predictions  require  an  accurate  flow  model  in  which  the  strong 
viscous-inviscid  coupling  is  modelled  correctly.  Boundary  layer  theory  correctly 
models  this  coupling  only  when  the  boundary  layer  is  closely  attached  to  the  body 
surface.  In  the  case  of  vortex  roll-up  on  the  lee  side  of  an  ogival  body  at  angle 
of  attack,  boundary  layer  theory  is  no  longer  an  adequate  model.  Although  a 
numerical  solution  to  the  full  Navier-Stokes  equations  would  provide  a correct  model, 
three-dimensional  Navier-Stokes  solutions  should  be  used  only  if  no  suitable 
alternative  exists.  A suitable  alternative  should  lead  to  accurate  solutions,  but 
should  not  be  limited  by  the  large  running  time  and  storage  requirements  associated 
with  three-dimensional  Navier-Stokes  solutions. 

The  flow  over  an  ogival  body  at  incidence  in  a supersonic  free  stream  can  be 
decomposed  into  two  distinct  parts.  First,  there  is  a streamwise  part  which  is 
primarily  determined  by  inviscid  mechanisms.  Then  there  is  a second  part  which  is 
orthogonal  to  the  first.  This  latter  part  accounts  for  the  secondary  flow  around 
the  body  and  tends  to  be  viscous  dominated.  Experimental  evidence  supports  the 
above  decomposition  (see  Refs.  1,  2,  3 and  l).  At  zero  angle  of  attack  a symmetric 
flow  field  develops  about  the  body,  in  which  no  cross  flows  exist.  As  the  angle 
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of  attack  is  increased  from  zero,  a cross  flow  pattern  begins  to  develop.  As  the 
angle  is  further  increased,  a cross  flow  adverse  pressure  gradient  develops  and 
becomes  sufficiently  strong  to  cause  flow  separation  and  then  a pair  of  vortices 
symmetric  about  the  lee  symmetry  plane  as  depicted  in  Fig.  1.  These  vortices  can 
have  imbedded  shocks  associated  with  them.  The  flow  as  depicted  in  Fig.  1 has  been 
verified  for  both  laminar  and  turbulent  flows  for  speeds  in  the  supersonic  and  hyper- 
sonic regimes.  The  available  experimental  data  also  indicates  that  flow  separation 
in  the  axial  direction  usually  is  not  associated  with  the  vortex  development. 

The  majority  of  previous  attempts  to  predict  flows  of  this  type  are  numerical 
solutions  of  the  inviscid  flow  equations,  the  three-dimensional  boundary  layer 
equations  or  both.  Furthermore,  the  geometry  was  invariably  restricted  to  conical 
axially  symmetric  bodies.  Such  analyses  cannot  predict  the  vortices  or  be  applied 
to  problems  with  a nontrivial  body  geometry. 

The  most  common  procedure  is  a numerical  solution  of  the  time-dependent 
inviscid  flow  equations  to  obtain  a steady  state  flow  field  (e.g..  Refs.  5 and  6). 
Other  available  inviscid  techniques  include  the  inverse  method  (Ref.  7)  and  the 
method  of  integral  relations  (Ref.  8).  MacCormack  and  Warming  (Ref.  9)  have 
recently  surveyed  the  available  inviscid  computational  procedures.  Viscous  forces 
may  in  part  be  accounted  for  by  making  use  of  the  inviscid  pressure  distribution  to 
solve  the  three-dimensional  boundary  layer  equations  (see  Ref.  10).  The  range  of 
applicability,  however,  places  an  unacceptable  restriction  on  the  angle  of  attack. 

At  significant  angles  of  attack,  there  is  a need  to  model  the  complex  flow 
field  about  the  body  with  a minimal  amount  of  effort.  This  effort  is  complicated 
by  difficulties  associated  with  the  synthesis  of  inviscid  flow  analysis  and  boundary 
layer  theory  into  a cohesive  method  of  analysis.  Among  these  difficulties  are  the 
lack  of  applicability  of  three-dimens inal  boundary  layer  theory,  a means  for 
patching  or  interfacing  boundary  layer  and  rotational  inviscid  flow  regions,  and 
the  treatment  of  interaction  between  viscous  and  inviscid  flow  regions. 

In  efforts  to  circumvent  the  above  difficulties,  Patankar  & Spalding  (Ref.  11), 
Caretto,  Curr,  & Spalding  (Ref.  12),  and  Briley  (Ref.  13)  constructed  numerical 
methods  for  solving  approximate  governing  equations  which  can  be  viewed  as  a general- 
ization of  three-dimensional  boundary  layer  theory.  In  these  studies,  solutions  were 
computed  for  laminar  incompressible  flow  in  straight  ducts  with  rectangular  cross 
sections.  The  governing  equations  were  solved  by  integrating  in  a primary  flow 
coordinate  direction  while  retaining  viscous  stresses  in  both  transverse  coordinate 
directions  as  opposed  to  only  one  direction  for  three-dimensional  boundary  layer 
theory.  In  addition,  the  pressure  gradient  terms  were  given  special  treatment 
to  permit  solution  by  forward  marching  integration.  Subsequently,  this  general 
approach  has  been  used  to  compute  laminar  incompressible  flow  in  helical  tubes  by 
Patankar,  Pratap  & Spalding  (Ref.  lU).  A predictor/corrector  solution  procedure 
has  been  developed  Lin  and  Rubin  (Refs.  15  and  16)  to  solve  the  approximate  three- 
dimensional  compressible  Navier-Stokes  equations.  The  numerical  technique  is 
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implicit  in  one  transverse  direction  and  iterative  in  the  other.  Hellivell  and 
Lubard  (Ref.  17),  Rakich  and  Lubard  (Ref.  18),  and  Line  and  Rubin  (Ref.  l6)  have 
applied  this  method  to  flows  over  both  sharp  and  spherically  blunted  cones  at  angle 
of  attack. 

In  companion  studies,  Briley  & McDonald  (Ref.  19)  and  McDonald  & Briley  (Ref.  20) 
have  developed  implicit  techniques  for  systems  of  coupled  nonlinear  equations.  These 
techniques  were  applied  by  McDonald  & Briley  (Ref.  20)  to  the  computation  of  laminar 
supersonic  flow  in  rectangular  jets.  Subsequently,  the  laminar  incompressible 
straight-duct  analysis  of  Briley  (Ref.  13)  and  the  improved  numerical  techniques 
of  McDonald  & Briley  (Ref.  20)  for  compressible  flows  were  extended  by  Briley  & 
McDonald  (Ref.  2l)  and  Eiseman,  McDonald,  and  Briley  (Ref.  22)  into  a method  for 
computing  subsonic  turbulent  flow  in  curved  ducts.  The  present  study  represents 
a further  generalization  of  the  latter  method;  it  encompasses  general  coordinates, 
mixed  subsonic  and  supersonic  flow,  and  shock  waves;  it  is,  in  fact,  a continuation 
of  the  study  initiated  by  Eiseman  and  Levy  (Ref.  23). 

The  basic  geometry  in  the  ogival  body  problem  is  determined  by  both  the  body 
itself  and  its  bow  shock.  Unlike  the  body  shape,  the  bow  shock  is  not  known  in 
advance,  but  must  be  determined  as  part  of  the  solution.  Since  the  region  nearest 
the  shock  is  dominated  by  convective  forces,  and  the  shock  is  treated  as  a discon- 
tinuity, it  is  quite  sufficient  to  perform  a local  inviscid  analysis  in  that  region; 
and,  thereby  to  determine  the  shock  location  one  step  in  advance  of  the  fully 
viscous  solution.  The  shock  location  is  calculated  numerically  in  terms  of  local 
extensions  of  the  existing  coordinate  system.  Next  the  shock  surface  is  extended 
in  a sufficiently  smooth  manner  so  that  suitable  coordinates  can  be  generated  for 
the  advancement  of  the  fully  viscous  solution.  For  this  purpose,  a least  squares 
surface  fitting  procedure  is  used. 

A major  part  of  the  solution  procedure  is  the  viscous  approximation  to  the 
Navier-Stokes  equations.  Specifically,  the  tensor  form  of  the  Navier-Stokes 
equations  is  approximated  in  a coordinate  independent  manner  to  produce  an  initial 
value  problem.  The  approximation  is  obtained  from  a neglect  of  viscous  stresses, 
and  hence  diffusion,  in  an  assumed  primary  flow  direction.  The  approximation  can 
be  viewed  as  a generalization  of  the  previous  work.  The  primary  flow  direction  is 
assumed  to  be  given  by  some  smooth  vector  field.  Since  the  specification  of  any 
vector  field  is  independent  of  coordinates,  the  approximation  of  the  tensor  form 
of  the  Navier-Stokes  equations  is  also  independent  of  coordinates.  As  a matter  of 
convenience,  however,  the  primary  vector  field  is  often  chosen  to  be  the  vector  field 
associated  with  a given  coordinate  direction  of  a given  system  of  coordinates. 

The  final  phase  in  the  development  process  of  a general  code  is  the  application 
and  extension  to  problems  of  ever  increasing  difficulty.  This  phase  can  only  be 
entered  once  the  underlying  structure  of  the  general  code  has  been  established  and 
associated  coding  errors  have  been  removed.  In  short,  a general  code  must  be  brought 
up  to  a reasonable  level  of  development  before  it  can  be  applied  and  extended  to 
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cover  a sequence  of  problems  leading  to  a particular  objective.  In  the  case  of  PEPSIG 
the  desired  level  of  development  was  recently  achieved.  Applications  and  extension 
of  the  code  have  now  been  made  for  both  internal  and  external  flow  problems.  A 
sequence  of  internal  flow  problems  has  lead  to  successful  calculations  of  subsonic 
flows  in  ducts  with  superelliptically  varying  cross  sections  under  a contract  with 
NASA  Lewis  Research  Center.  The  objective  for  the  external  flow  problems  under  the 
present  effort  was  the  computation  of  the  supersonic  flow  over  an  ogival  body  at 
angle  of  attack.  This  objective  involved  the  solution  of  mixed  supersonic  and  sub- 
sonic flow;  specifically,  a subsonic  layer  is  bounded  by  a solid  object  and  a super- 
sonic flow. 

Before  significant  progress  can  be  made  on  the  overall  program  goal,  the  code 
must  be  able  to  accurately  model  shock  waves,  subsonic  boundary  layers,  and  super- 
sonic flow  with  no-slip  boundary  conditions.  Both  shock  waves  and  subsonic  boundary 
layers  have  been  modelled  adequately  and  the  results  will  be  displayed  along  with 
preliminary  results  for  supersonic  flow  with  no-slip  boundary  conditions.  This 
last  case  involving  mixed  subsonic  and  supersonic  flow,  however,  has  proven  to  be 
quite  difficult  to  solve  to  our  desired  level  of  accuracy.  Most  investigators 
(Refs.  16-18)  have  noted  the  development  of  "departure  solutions"  as  an  approximate 
set  of  equations  governing  this  mixed  flow  is  marched  in  the  axial  direction. 

Suggested  methods  to  eliminate  these  unwanted  departure  solutions  have  been  to  open 
up  the  mesh  spacing  to  meet  a certain  lower  bound,  to  specify  a pressure  gradient  in 
the  axial  direction,  and  to  carefully  prescribe  an  initial  profile.  Although, 
a lower  bound  on  the  axial  mesh  may  be  required  to  obtain  a stable  marching 
solution  of  certain  stiff  differential  equations  (Ref,  24),  it.  is  disturbing  when 
accuracy  is  also  desired  as  is  the  case  that  is  envi  > toned  for  cK.  ogival  body’ 
with  curvature.  In  any  case,  when  the  axial  mesh  ia  ru.  .'aseu  flow  over  a 
wedge,  in  the  present  cases  only  a decrease  in  th<  v».  of  growth  of  departure 
solutions  is  observed.  The  symptoms  of  departs  , Monr  however,  developed 
immediately  as  indicated  by  the  development  of  *.<.■ . ,'se  ve:  ^ure  gradients  in  the 
immediate  vicinity  of  the  wall.  It  is  also  curious.  note  that,  up  to  a certain 
point,  the  present  solutions  had  many  qualitatively  correct  features  when  viewed 
on  isolated  cross  sections.  This  predicament  tends  to  make  one  concerned  that  the 
previously  developed  remedies  to  suppress  departure  solutions  are  in  fact,  case 
dependent;  it  also  points  to  a need  for  an  additional  general  development.  Recently, 
the  problem  of  departure  solutions  has  been  more  thoroughly  investigated  in  Refs. 

25  and  26.  These  investigators  observed  that  supersonic  pressure  waves  would 
propagate  above  the  subsonic  portion  of  the  boundary  layer  and  influence  it  in  the 
upstream  direction  causing  a pressure  rise  which,  if  left  alone,  would  cause  the 
solution  to  depart  from  the  underlying  physics.  Their  remedy  was  to  specify  an 
additional  condition  at  the  sonic  line;  namely,  that  the  velocity  is  parallel  to 
the  wall  angle.  While  this  though  proved  to  be  more  successful  than  that  of  previous 
investigators,  when  applied  in  the  present  instance,  entirely  satisfactory  results 
were  not  obtained.  In  addition,  two  variants  of  this  technique  were  tried  in  the 
present  effort.  Firstly,  an  inviscid  velocity  angle  was  specified  at  the  sonic  line 
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and  secondly,  an  extrapolated  angle  was  tried.  The  latter  variant  resulted  in  a 
slight  improvement . Additional  developments  centered  around  new  specifications  of 
the  axial  pressure  gradient  in  the  subsonic  region  and  on  the  treatment  at  the  sonic 
line.  To  date  the  results,  while  improving,  are  still  not  satisfactory.  Further 
work  on  this  particular  problem  is  obviously  required. 
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THE  GENERAL  FORMATION  OF  AN  INITIAL  VALUE  PROBLEM  FOR  STRONGLY  CONVECTIVE  FLOWS 


Central  to  the  present  analysis  is  the  formulation  of  approximate  governing 
equations  which  can  be  solved  by  forward  marching  integration  in  the  direction  of 
a 'primary  flow".  The  entire  flow  field  can  thereby  be  obtained  by  a sequence  of 
essentially  two-dimensional  calculations.  This  feature  of  the  method  results  in 
a substantial  saving  of  computer  time  and  storage  compared  to  that  which  would  be 
required  for  solution  of  the  full  Navier-Stokes  equations.  The  equations  are 
derived  in  a coordinate  independent  manner.  The  merits  of  such  a derivation  is  given 
Ref.  23.  A vector  field  that  reasonably  approximates  the  primary  flow  direction 
is  chosen  and  then  used  as  the  basis  for  an  approximation  of  the  stress  tensor. 

The  time-averaged  equations  are  written  in  general  conservation  law  form,  and  then 
the  approximate  stress  tensor  is  inserted  to  obtain  the  approximate  equations. 

Note  that  this  process  depends  only  on  the  choice  of  a primary  vector  field,  and 
not  on  the  particular  coordinate  system  used  for  the  numerical  solution.  The 
primary  vector  field  used  here  consists  of  the  tangent  vectors  to  a certain  family 
of  coordinate  curves  that  are  roughly  aligned  with  the  flow  geometry. 

The  governing  equations  are  derived  from  the  Navier-Stokes  equations  for  com- 
pressible flow  of  a viscous,  perfect  gas.  In  conservation  law  form  (Refs.  27-28) 
and,  in  general  curvilinear  coordinates  (y1,  y^,  y3),  the  continuity  equation  is 
given  by 


(„‘Vg)  = o 
dt  ' ' 


the  momentum  equations , by 


If  constant  total  temperature  is  assumed,  and  energy  equation  is  not  required  and 
can  be  replaced  by 

pv[»‘«v"'l  (1C) 

with  the  added  assumption  of  a perfect  gas.  Otherwise,  an  energy  equation 

+»rir'ivr|v/5'].0  da] 
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oust  be  used.  In  the  above  (x1,  x2,  x^)  represents  fixed  cartesian  coordinates 
p,  density;  E,  energy;  k,  thermal. conductivity;  A and  B,  constants;  Vs  » 


velocity; 


9 * dettyj) « 


W 


m 


the  metrical  determinant;  and  t*J,  the  components 


of  the  stress  tensor  in  the  basis  5^  ^ei,ias  ne,tric,  the  components 

of  the  stress  tensor  are  given  by 


where 


r'i  s g'i  p + ak'i  vk  + 

dy 


(2a) 


and 


(2b) 


(2c) 


for  viscosity/*,  inverse  metric  g*^  , Kronecker  deltas  8lj  = 8,*  = 8j.  * and 
Chris toff el  symbols  * 


rk  = -i!!  I foim  ag^m  -*9jj) 
'i  2 ( dy  * dy 1 dym j 


(2d) 


In  all  of  the  above,  the  Einstein  summation  convention  is  assumed.  That  is,  match- 
ing upper  and  lower  indices  are  to  be  summed  from  1 to  3 unless  otherwise  stated. 


It  is  assumed  that  for  high  Reynolds  number,  viscous  effects  are  negligible 
except  in  thin  layers  near  the  wadis,  and  thus  boundary  layer  concepts  can  be 
enployed  to  examine  relative  importance  of  viscous  terms  in  the  governing  equations. 
Consequently  viscous  terms  which  are  considered  important  for  boundary  layer  flow 
on  walls  are  retained;  other  viscous  terms  are  neglected.  In  this  sense,  the  pre- 
sent approach  can  be  regarded  as  a natural  extension  of  three-dimensional  boundary 
layer  theory.  Unlike  conventional  boundary  layer  theory,  however,  the  approximate 
equations  are  to  be  applicable  in  the  inviscid  flow  region  as  well  as  the  viscous 
region  and,  thus,  no  approximations  are  made  for  inviscid  terms  other  than  those 
to  be  used  for  the  pressure  field  in  subsonic  flow.  A detailed  account  of  the 
viscous  approximation  is  given  in  Ref.  23. 
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METHOD  OF  SOLUTION 


Background 

The  governing  equations  can  he  solved  following  the  general  approach  developed 
in  Ref.  23.  The  method  used  is  based  on  an  implicit  scheme  which  is  potentially 
stable  for  large  step  sizes.  Thus,  as  a practical  matter,  stability  restrictions 
which  limit  the  axial  step  size  relative  to  the  transverse  mesh  spacing  and  which 
become  prohibitive  for  even  locally  refined  meshes  (e.g.,  in  laminar  sublayers)  are 
not  a factor  in  making  the  calculations.  The  general  approach  is  to  employ  an  implicit 
difference  formulation  and  to  linearize  the  implicit  equations  by  expansion  about  the 
solution  at  the  most  recent  axial  location.  Terms  in  the  difference  equations  are  then 
grouped  by  coordinate  direction  and  one  of  the  available  alternating-direction  implicity 
(ADI)  or  splitting  techniques  is  used  to  reduce  the  multidimensional  difference  equations 
to  a sequence  of  one-dimensional  equations.  These  linear  one-dimensional  difference 
equations  can  be  written  in  bio ck-tri diagonal  or  a closely  related  matrix  form  and 
solved  efficiently  and  without  iteration  by  standard  block  elimination  techniques. 

The  general  solution  procedure  is  quite  flexible  in  matters  of  detail  such  as  the  type 
and  order  accuracy  of  the  difference  approximations  and  the  particular  scheme  for 
splitting  multidimensional  difference  approximations.  Based  on  previous  experience 
of  the  authors,  however,  it  is  believed  that  the  consistent  uBe  of  a formal  lineari- 
zation procedure,  which  incidentally  requires  the  solution  of  coupled  difference 
equations  in  most  instances,  is  a major  factor  in  realizing  the  potential  favorable 
stability  properties  generally  attributed  to  implicit  difference  schemes. 


Numerical  Techniques 

As  an  outline  of  the  particulars  of  the  numerical  method,  the  treatment  of  the 
continuity  equation  is  considered,  as  this  is  the  simplest  equation,  and  yet  this 
discussion  will  cover  most  aspects  of  the  method.  The  flow  region  is  discretized  by 
grid  points  having  spacings  Ax  * Ay3,  Ay  = Ayl,  and  Az  * Ay2  respectively.  Pro- 
visions for  nonuniform  grid  spacing  can  be  introduced  as  needed.  The  subscripts 
i,  J and  superscripts  n are  grid  point  indices  associated  with  yey^,  z=y^,  and 
x = y3,  respectively.  Thus,  . denotes  $(xn,yi,z.)  where  4>  can  represent  any  of 
the  dependent  variables.  The  subscripts nare  frequently  omitted  if  clarity  is 
preserved,  so  that  is  equivalent  to  j,.  For  convenience,  the  following  shorthand 
difference  operator  notation  is  used  for  derivative  difference  formulas: 


,-2+M+*r+u]/w' 


(3a) 


(3b) 
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with  analogous  definitions  for  ( 8* . Here , a parameter  a has  "been  introduced 
(0  set  SI)  so  as  to  permit  continuous  variation  from  backward  to  forward  differences. 
The  standard  central  difference  formula  is  recovered  for  a s 1/2  . Throughout  the 
following  discussion,  it  is  assumed  that  the  solution  is  known  at  Xn  and  is  desired 

at*"  + l. 

Consider  the  laminar  continuity  equation  with  v5  *U,VI*V,  v* 


•£-  ‘ (pu J)  + -i-  (pvJ)+  ~ (pwJ)  = 0 
ox  oy  oz 


(M 


where  J • Equation  (U)  is  differenced  in  the  x or  marching  direction  as 
follows : 


where 


^n  + /9_  ♦n+'+/9$n 

* ~ FT — 


ond  0- 


0 

l+/3 


Here,  a parameter  8 has  been  introduced  bo  as  to  permit  a variable  centering  of  the 
scheme  in  the  x direction.  Equation  (5)  produces  a backward  difference  formulation  for 
8*0  and  a Crank-Nicolson  formulation  for  8*1.  The  dependent  variables  in  Eq.  (5) 
are  linearized  in  the  manner  of  Refs.  (19)  and  (20)  by  expansion  about  the  solution 
at  at0.  Here,  a first-order  accurate  linearization  is  used  for  the  x derivative,  and  the 
result  is 


9 


RT7-912536-8 


. 


| 


» 

i 


r ^nun  + 1 + pti  + 1 u".2^nun 

r jn+,+^Jn' 

l Ax 

l 

(6) 


+ 77a{j"  + '»y  VvB+,V+l*n-^*"l  + (li {•)"*Vnyn+'V  + lvn-/'v") 


/3 


j"  + '*,(,"w"+,+/+la,W)  + (|i)  (A>"+,VV-A«n) 
+ p [js,  <a>">  +(|f)V«-"]}  *0 


|J 

| 


. 


j 

j 

1 


I 


Neglecting  for  the  moment  the  turbulent  Reynolds  stresses,  and  eliminating  the 
pressure  gradients  via  the  equation  of  state,  the  procedure  outlined  above  for  the 
continuity  equation  can  be  employed  to  derive  linear  Implicit  difference  approxi- 
mations analogous  to  Eq.  (6)  for  the  three  momentum  equations.  Furthermore,  the 
resulting  difference  approximations  can  be  grouped  by  coordinate  direction  and 
written  in  the  following  compact  linear  matrix  difference  operator  notation  as 

A /*n  +1  n\  An+I  . _n+l  „ 

^ | « Dy1  + + S (T) 

where  $ is  a column  vector  containing  the  dependent  variables  p,  u,  v,  u,  and  A 
is  a square  (4x1*)  matrix.  Dy  and  D are  4x4  matrices  containing  elements  which  are 
themselves  spatial  difference  operators  for  the  y and  z directions,  respectively.  S 
is  a column  vector  reserved  for  any  source  terms  which  may  be  present.  The  matrices 
A,  D , and  Dz  contain  only  quantities  which  are  known  from  a computational  viewpoint. 
Equation  (7)  is  linear  in  4>n+1. 
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Solution  of  the  Split  Difference  Equations 


The  coupled  set  of  linear  implicit  difference  equations  arising  along  rows  of 
grid  points  during  each  step  of  the  ADI  solution  procedure,  together  with  the 
prescribed  boundary  conditions,  can  be  written  in  a form  having  the  matrix  structure 
given  by  a modified  block  tridiagonal  form.  Neglecting,  for  the  moment,  the 
modifications  that  resulted  from  boundary  conditions,  a direct  solution  by  standard 
block  elimination  techniques  (cf. , Isaacson  & Keller,  Ref.  3l)  is  both  straightforward 
and  efficient.  The  precise  scheme  used  here  consisted  of  Gaussian  elimination  for  a 
simple  tridiagonal  system  (sometimes  called  the  Thomas  algorithm)  but  with  elements 
of  the  tridiagonal  matrix  treated  as  square  submatrices  rather  than  as  simple  coeffi- 
cients. The  required  inverses  of  diagonal  submatrices  were  obtained  by  a Gauss- 
Jordan  reduction.  The  additional  operations  necessary  to  include  the  nonblock- 
tridiagonal  elements  are  easily  incorporated  provided  the  original  block  tridiagonal 
coding  is  carefully  organized. 
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THE  COORDINATE  SYSTEM 

J 

The  Construction  of  Tube-like  Coordinates 

3 

Tube-like  coordinates  will  be  constructed  in  general  to  provide  a natural 
setting  for  the  study  of  flows  within,  between,  or  outside  of  a set  of  prescribed 
tubes.  The  prescribed  boundary  tubes  then  become  coordinate  surfaces,  and,  as  a 
result,  the  specification  of  fluid  dynamic  boundary  conditions  is  greatly  simplified. 
Although  the  equations  of  motion  contain  more  terms  than  for  a cartesian  system, 
this  does  not  add  excessively  to  the  run  time  of  a program.  In  addition,  there 
must  be  some  control  over  the  resolution  of  regions  near  bounding  tubes  so  that  the 
effects  of  wall  curvature  and  the  growth  of  attached  boundary  layers  can  be  ade- 
quately treated.  Such  controls  are  obtained  from  the  specification  of  coordinate 
distribution  functions  which  shall  appear  only  as  parameters  in  the  basic  geometric 
construction  of  the  coordinates.  The  basic  geometry  of  the  bounding  tubes  then 
provides  the  intrinsic  constraints  upon  the  coordinate  construction.  Since  the 
primary  goal  is  the  computation  of  fluid  flows  within  nontrivial  geometries  and  not 
the  development  of  coordinate  systems  per  se,  the  coordinates  will  be  kept  as  simple 
as  possible,  given  the  desired  generality. 


Considering  various  past  successes  of  two-dimensional  conformal  mappings  to 
obtain  coordinates,  one  might  naturally  wish  to  obtain  similar  transformations  for 
three-dimensions.  Unfortunately,  there  is  no  three-dimensional  theory  of  conformal 
transformations  analogous  to  complex  variables,  and  consequently,  in  three  dimen- 
sions one  is  left  with  a complicated  system  of  partial  differential  equations  which 
generally  would  require  numerical  solution.  To  circumvent  the  considerable  computa- 
tional labor  required  for  solution  of  such  equations,  a constructive  process  is 
used  for  the  development  of  tube-like  coordinates. 

The  first  step  in  the  construction  of  tube-like  coordinates  is  to  create  a 
suitable  family  of  two-dimensional  surfaces  which,  in  some  sense,  are  transverse  to 
a given  centerline.  If  orthogonal  coordinates  are  desired,  then  these  surfaces  would 
have  to  bend  and  flex  as  the  tube  would  undergo  changes  in  cross  section  at  differ- 
ent centerline  positions,  m addition  to  the  problem  of  constructing  transverse 
surfaces  which  bend  and  flex,  there  is  also  the  problem  of  constructing  an  orthogonal 
grid  on  a surface  which  has  variations  in  Gaussian  curvature,  and  hence,  is  not  flat. 
This  second  problem,  in  fact,  requires  a more  complicated  construction  than  the 
first  which  in  itself  is  not  easy.  Thus,  the  sheer  magnitude  of  the  work  involved 
in  the  construction  of  orthogonal  coordinates  certainly  would  remove  the  desire 
for  their  use  in  fluid  dynamic  problems  which  undoubtedly  would  require  less  com- 
putation in  nonorthogonal  coordinates  than  in  the  construction  of  an  orthogonal 
system  alone.  By  contrast,  if  the  transverse  surfaces  are  selected  to  be  two-dimen- 
sional planes,  then  the  construction  of  coordinates  is  greatly  simplified  while  the 
fluid  dynamic  computation  is  only  marginally  different  due  to  coordinate  nonortho- 
gonality. Consequently  the  coordinate  system  that  we  shall  construct  will  have 
planar  transverse  surfaces . 
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Since  each  planar  transverse  surface  is  a linear  subspace  of  the  real 
three-dimensional  Euclidian  vector  space  r3,  any  such  plane  can  be  completely 
specified  by  any  two  spanning  linearly  independent  vectors  in  R^.  The  specifica- 
tion of  the  planar  family  of  transverse  surfaces  is  then  a result  of  a construction 
of  two  vector  fields  along  a given  centerline  curve  in  R . The  origin  of  each  plane 
is  chosen  to  coincide  with  the  associated  centerline  point.  (See  Fig.  2.)  To 
assure  that  the  planes  are  always  indeed  transverse,  it  will  be  assumed  that  they 
are  orthogonal  to  the  centerline  at  their  origins.  Ultimately,  tube-like  surfaces 
will  be  generated  by  loops  about  the  planar  origins  which  deform  in  some  way  as  we 
move  along  the  centerline  curve;  in  general,  these  tube-like  surfaces  will  not 
intersect  the  transverse  planes  orthogonally.  Thus,  only  the  centerline  direction 
determines  the  transverse  nature  of  the  cross  sectional  planes.  Specifically, 
the  centerline  tangent  vectors  form  a vector  field  which,  at  each  point,  is  ortho- 
gonal to  the  plane  of  the  two  transverse  vectors,  and  thus  each  centerline  point 
carries  a triple  of  linearly  independent  vectors.  By  the  Gram-Schmidt  orthogonali- 
zation  procedure,  each  such  triple  of  vectors  can  be  made  into  an  orthogonal  set, 
and  hence,  an  orthonormal  set  which  is  simply  called  a frame.  Thus,  tube-like 
coordinate  systems  are  constructed  from  a specified  centerline  curve  and  an  asso- 
ciated frame  field.  Row  the  basic  question  is  whether  there  is  a canonical  con- 
struction of  tube-like  coordinate  systems  from  either  a given  centerline  or  a given 
frame  field.  From  the  theory  of  space  curves  (Ref.  32),  it  is  well  known  that  for 
positive  curvature  and  specified  torsion  there  is  a local  one-to-one  correspondence 
between  frame  fields  and  space  curves  which  pass  through  a given  point.  Thus,  for 
nonzero  curvature , the  centerline  space  curve  has  a canonical  frame  field  which  is 
known  as  the  Frenet  frame.  Consequently,  the  coordinates  will  be  derived  from  the 
Frenet  frame  when  it  exists.  At  centerline  points  of  zero  curvature,  the  Frenet 
frame  is  degenerate  and  must  be  treated  specially. 

Once  the  Frenet  frame  of  the  space  curve  has  been  established,  the  unit 
normal  and  binormal  vectors  Vg  and  at  each  point  of  determine  a transverse 
plane  orthogonal  to  the  unit  tangent  vector  (See  Fig.  :3. ) Relative  to  any  such 

transverse  plane,  these  vectors  are  also  the  standard  orthonormal  basis.  Conse- 
quently, we  can  examine  the  plane  separately  from  the  curve,  V,  which  will  only 
appear  as  the  point  at  the  origin.  In  two  dimensional  functional  terminology,  the 
unit  normal  direction  can  be  considered  as  the  abscissa  and  the  unit  binormal  as  the 
ordinate;  or  more  simply,  as  x and  y axes,  respectively.  Since  the  tube-like 
coordinates  are  to  be  generated  from  some  family  of  tubes  encasing  the  space  curve, 
Y,  a cross-sectional  cut  by  a transverse  plane  produces  within  the  plane  a family 
of  loops  about  the  origin.  We  shall  assume  that  each  loop  is  representable  by  a 
strictly  monotone  radial  function  of  angle.  In  this  regard,  a polar  type  of  descrip' 
tion  is  the  most  suitable.  But,  of  course,  the  loops  are  usually  more  complicated 
than  circles,  and  thus,  we  must  replace  the  radius  by  a function  L of  both  radial 
and  angular  variables  r and  e.  Furthermore,  when  noncircular  loops  bound  a cross 
section  of  fluid,  there  are  regions  of  varying  wall  curvature.  In  a numerical 
solution,  it  is  desirable  to  put  proportionately  more  mesh  points  in  regions  of 
higher  curvature  than  in  regions  of  less  curvature . Consequently,  an  angular 
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distribution  function,  6,  is  a good  replacement  for  the  simple  angular  specification, 
8,  of  simple  polar  coordinates.  The  net  result  is  a generalization  where  polar 
coordinates  are  replaced  by  a pseudo-radius,  r and  pseudo-angle,  8.  Since  the  loops 
generally  vary  from  transverse  plane  to  transverse  plane,  the  pseudo-radii  and  angles 
must  also  be  functions  of  axial  location,  t,  on  the  centerline  space  curve,  "ytt). 

Since  the  normal  and  binormal  directions  are  usually  functions  only  of  the  centerline 
curve,  y(t),  our  loops  may  have  symmetries  that  do  not  reflect  about  either  of  these 
Frenet  directions.  Since  the  use  of  known  symmetries  is  a great  simplification  in 
most  problems,  we  need  an  option  which  allows  one  to  define  axes  that  can  be  aligned 
in  an  optimal  way.  This  option  is  easily  established  from  the  specification  of  a 
function,  fl(t),  which  is  a rigid  rotation  relative  to  the  normal -binormal  directions. 
To  bring  this  development  of  tube-like  coordinates  within  the  framework  of  the 
preceding  tensor  derivations,  we  shall  use  the  notation,  y^=«8,  y^=r  and  y3*t  for 
pseudo- angular,  pseudo-radial,  and  axial  variables.  In  this  notation,  we  have  thus 
far  developed  (1)  a length  factor,  L=L(y1,  y2,  y3),  which  is  a generalization  of 
radius,  (2)  an  angular  distribution  function,  JMjjy1,  y3)  which  is  a generalization 
of  angle,  (3)  a rotation  function,  (tn(y3),  and  (4)  the  Frenet  frame,  (Vj,  V,,,  Vl)= 
(Vi(y3),  T2(y3),  V^y^))  upon  which  the  coordinates  are  built.  That  the  length 
factor,  L,  and  the  angular  distribution  function,  0,  give  us  a generalization  of 
polar  coordinates  is  obvious  since  polar  coordinates  are  easily  retrieved  by  taking 
L(y*,  y2,  y3)  = y2  and  Sty1,  y3)  * y*.  It  is  also  worth  noting  that  the  angular 
distribution  function,  0,  was  chosen  to  be  independent  of  pseudo-radius,  y2. 

Although  it  is  not  immediately  evident,  we  have  removed  a considerable  amount  of 
potential,  computational  complexity  in  the  process  of  obtaining  metric  information 
by  limiting  the  number  of  derivatives  which  must  be  computed.  Furthermore,  there 
is  no  real  loss  of  flexibility  in  the  construction  of  angular  distribution  functions. 
Since  most  commonly  used  analytic  descriptions  of  loops  are,  in  fact,  controlled  by 
a collection  of  parameters  which  depend  only  on  axled  location,  y3,  a knowledge  of 
only  these  parameters  is  often  sufficient  for  the  construction  of  the  angular  dis- 
tribution function.  For  example,  if  the  loops  were  to  consist  of  a family  of  con- 
centric homogeneous  ellipses,  then  the  major  and  minor  axes  of  the  outermost  ellipse 
would  form  a collection  of  two  such  parameters. 

With  the  above  functions  and  the  Frenet  frame,  the  class  of  tube-like  coordinates 
comes  directly  out  of  the  transformation 

7=7+ cos<#>+73  sin$} 

which  transforms  curvilinear  coordinates,  y « (j4,  y2,  y3)  into  cartesian  coordinates 
¥ • (x1,  x2,  x3)  where 

<My,.ys)*®(y,iy5)*  fi(  y5) 
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At  each  transverse  location,  y3,  the  space  curve  vector,  ^ translates  the  origin 
to  the  space  curve.  At  a given  pseudo-angle,  y^-,  a unit  vector,  V2  cos  cp  + sin  <P> 
determined  by  the  sum,  (p  = 0 + fi  of  the  radial  distribution  function,  e,  and  the 
transverse  rotation,  0.  This  unit  vector  sweeps  out  a full  3&0  degs  in  the  trans- 
verse plane  as  y1  passes  through  all  of  its  values.  Hence,  we  could  call  this  a 
direction  pointer  for  the  transverse  plane.  When  this  direction  pointer  is  scaled 
by  the  length  factor,  L,  we  obtain  a point  of  our  transformation.  Since  the  length 
factor  depends  on  all  three  variables,  any  set  of  tube-like  surfaces  can  be  obtained 
provided,  of  course,  that  loops  are  representable  by  a strictly  monotone  radial 
function  of  angle  and  also  that  no  two  transverse  cross  sections  are  allowed  to 
intersect. 

In  a geometric  setting,  the  transformation  is  really  an  embedding  of  tube-like 
coordinate  systems  into  three  dimensional  Euclidian  space.  An  illustration  is 
provided  in  Fig.  3 • From  the  transformation,  it  is  also  easy  to  see  that  the  surfaces 
of  constant  y3  are  the  transverse  planes,  the  surfaces  of  constant  pseudo-angle,  y^, 
are  ruled  surfaces  generated  from  the  centerline  curve,  "y,  and  the  surfaces  of  con- 
stant pseudo-radius,  y^,  are  just  the  concentric  tubes  about  the  space  curve,  ‘y. 

Separate  illustrations  of  these  various  coordinate  surfaces  are  given  in  Figs.  4a, 

4b,  and  4c , respectively.  The  metric  data  for  this  transformation  is  given  in  Ref.  23.  j 
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The  Length  Factor 

Within  the  structure  of  tube-like  coordinates,  the  length  factor  contains  the  in- 
formation needed  both  for  the  specification  of  basic  geometry  and  for  the  distributional 
control  of  the  flow  region.  The  distributional  control  can  be  easily  implemented  by 
the  use  of  pseudo-radial  and  angular  distributions  as  merely  parameters  in  the  con- 
struction of  the  length  factor.  In  this  way,  the  basic  geometry  can  be  treated  sepa- 
rately from  the  distributional  aspects.  The  point  of  separation  becomes  especially 
evident  when  the  process  of  length  factor  construction  is  broken  down  into  a number 
of  stages.  Since  the  bounding  tubes  can  be  smoothly  generated  from  bounding  loops 
within  each  transverse  plane,  it  is  sufficient  to  temporarily  restrict  our  analysis 
to  a consideration  of  regions  between  bounding  loops  within  a transverse  plane.  For 
our  tube-like  coordinates,  we  have  implicity  assumed  that  the  nondegenerate  bounding 
loops  do  not  pass  through  the  origin  and  are  contractable  along  radial  lines  emanating 
from  the  origin.  Thus,  we  do  not  allow  bounding  loops  to  intersect  a radial  line 
more  than  once.  Consequently,  each  bounding  loop  can  be  expressed  in  terms  of  polar 
coordinates  (r,«)  as  the  product  of  a single  valued  radial  function  Fi(0)  and  the  unit 
vector  (cos  6,  sin  0).  With  this  polar  description  the  contraction  process  effectively 
reduces  to  a one-dimensional  process.  Specifically  it  is  a process  between  the  coef- 
ficients of  the  unit  vector  (cos  0,  sin  0).  For  a given  set  of  loops  Y^, ...,Yk  any 
sufficiently  smooth  interpolation  process  between  the  coefficients  will  be  satisfactory. 
If  we  assume  that  no  two  tubes  join  within  the  flow  region,  then  the  flow  region  is 
divisible  into  subregions  with  no  more  than  two  bounding  tubes.  The  interpolation  j 

process  for  two  bounding  loops  v1,  and  iB>  therefore,  all  that  is  usually  needed; 
and  is  given  by  the  simple  homotopy 
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H(0,r ) = L(fl,r ) (cos  0,  Sin®)  (11a) 

with  length  factor 

l_(^,r)=  rF,(©)  +(l-r)F2  (0)  (lib) 

which  takes  y2(e)  = H(0,o)  uniformly  and  smoothly  into  y^(e)  = H(e,l)  as  r goes  from 
0 to  1.  (See  Ref.  33.).  This  is  illustrated  in  Fig.  5.  if  the  loop  y2  is  degenerate 
then  the  coefficient  Fg(e)  vanishes  and  the  length  factor  reduces  to  L(e,r)  = rF.(B), 
we  thus  have  the  cross  section  of  a duct  generated  by  one  loop.  By  the  continuity 
of  L as  a function  of  y2,  bhe  duct  generated  by  one  loop  y^  can  be  considered  as  a 
limit  of  annular  type  regions  between  loops  and  y^  and  y2  closes  tightly  upon  the 
origin.  This  concept  is  often  quite  useful  since  the  origin  in  coordinates  generated 
from  one  loop  suffer  from  the  same  singularity  problem  that  occurs  with  simple  polar 
coordinates.  This  singularity  can  be  circumvented,  however,  by  using  an  auxiliary 
loop  y2  which  is  near  enough  to  the  origin  to  create  a good  approximation  to  the 
original  region.  To  preserve  overall  accuracy  in  a numerical  computation,  ||y2I|  * 
F2(0)  must  be  less  than  the  numerical  truncation  error.  In  fact,  the  well-defined 
limiting  process  would  lead  one  to  believe  that  there  would  be  no  problem  at  all  in 
taking  j|y2ll  arbitrarily  small.  But  if  ||y2l|  is  taken  within  the  region  of  machine 
roundoff  error,  then  the  singularity  problem  may  reappear  by  default.  Consequently, 
it  is  best  to  choose  | |y2| | to  be  much  less  than  truncation  errors  but  greater  than 
roundoff  errors. 


The  final  stage  of  length  factor  construction  is  accomplished  by  a replacement  of 
the  polar  coordinates  r and  0 by  radial  and  angular  distribution  functions  R(r,t)  and 
0(0 ,t)  for  axial  location  t.  Now  since  R and  © are  to  be  the  actual  polar  locations 
of  a loop  we  must  reinterpret  r and  0 as  pseudo-radial  and  pseudo-angular  locations  on 
the  same  loop.  Within  this  context  the  two-tube  length  factor  becomes 


L(0,r,t)  * R(r,t)F,  (•($,  t)  )+[l-R  (r,t )]  F2  (e(0,t) ) 

and  the  associated  unit  vector  becomes 

(cos(9(0,t)+  a( t)),  $intt>(0,t)+  A( t)) 


(12) 


(13) 


where  the  rotation  fl(t)  of  the  Frenet  frame  has  been  included  for  completeness. 
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The  Construction  of  Bounding  Tubes 

As  the  reader  has  Just  seen,  the  construction  of  tube-like  coordinate  systems 
relies  upon  the  existence  of  smooth  bounding  tubes.  When  such  tubes  exist  at  the 
outset  of  a problem,  the  generation  of  coordinates  is  a straightforward  application 
of  the  development  above.  However,  if  the  bounding  tubes  are  unknown  at  the  outset, 
then  they  must  be  constructed  in  a smooth  enough  fashion.  In  such  cases  one  is 
often  given  a discrete  specification  of  a sequence  of  bounding  loops  which  first 
must  be  fit  with  a smooth  curve  and  then  must  be  Joined  to  form  a smooth  surface. 

This  circumstance  can  often  arise  out  of  the  convenience  associated  with  the  discrete 
specification  of  a surface  by  means  of  successive  cross  sectional  cuts.  If  this 
data  were  to  be  given  all  in  advance  of  the  intended  use  of  the  coordinate  system 
as  a whole,  then  the  smoothed  cross  sectional  loops  could  be  effectively  Joined  by 
fitting  them  together  with  splines  which  can  have  interior  knots  corresponding  to 
interior  loops.  However,  it  is  often  the  case  that  discrete  loop  data  is  only 
generated  one  station  in  advance  of  the  use  of  the  coordinate  system.  This  occurs, 
for  example,  when  the  problem  is  to  solve  for  the  viscous  flow  field  outside  of  an 
ogival  body  when  the  flow  is  predominately  supersonic.  While  the  ogival  body  sur- 
face is  known  in  advance,  the  location  of  the  bow  shock  is  not.  Thus  one  considers 
the  ogival  body  surface  as  the  unknown  outer  tube  which  one  wishes  to  use  for  the 
generation  of  tube-like  coordinates  to  allow  for  the  efficient  computation  of  the 
fully  viscous  flow  field.  Since  the  flow  field  in  a neighborhood  of  the  bow  shock 
is  largely  inviscid,  an  inviscid  explicit  solution  is  performed  iteratively  to 
obtain  discretely  the  geometry  of  the  bow  shock  at  one  station  in  advance  of  our 
known  soltuion  and  coordinate  system.  One  is  now  left  with  a fully  developed  bow 
shock  surface  preceeded  by  a discrete  cross  sectional  loop  of  bow  shock  data.  The 
problem  is  to  smoothly  fit  the  loop  and  then  smoothly  Join  the  result  to  obtain  a 
smooth  extension  of  the  surface.  Since  fluctuations  may  arise  from  the  discrete 
generation  of  the  bow  shock  data,  a least  squares  spline  procedure  is  used  to  fit 
the  loop  with  smoothness  up  to  three  continuous  derivatives.  This  type  of  least- 
squares  procedure  has  the  distinct  advantage  of  accurately  representing  the  surface 
normal  curvature  along  the  loop.  Now  one  has  a loop  y (e)  at  level  n+1  and  a 
bounding  tube  ending  on  a loop  yn(e)  at  level  n where  S&ere  are  known  derivatives 
in  the  axial  direction  t.  One  then  attaches  a surface  which  extends  the  tube 
from  7n  to  yn+l  with  the  smoothness  of  three  continuous  derivatives.  The  extension 
is  accomplished  with  the  tensor  product  form 

p(0,t)sZ*'  fi  (t)  X_  . \(B)  {lk) 

j<0  1 

which  takes  information  back  to  loop  y p.  At  the  beginning  r must  be  1 since  there 
are  only  two  available  loops.  The  process  continues  with  r increasing  until  a 
desired  maximum  r value  is  obtained.  From  there  on  r is  assumed  to  be  constant. 
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Distribution  Functions 

When  partial  differential  equations  are  discretized  in  terms  of  differences, 
the  derivatives  are  replaced  in  some  fashion  by  difference  quotients.  A simplifi- 
cation then  leads  to  the  difference  equations  that  we  solve.  Implicitly  in  the 
discretization,  however,  is  the  assumption  that  derivatives  are  accurately  estimated 
by  secant  lines.  But  then  the  exact  solution  may  experience  drastic  variations  in 
a short  distance.  Such  solutions  are  said  to  have  large  gradients.  Ih  regions 
where  the  gradients  cure  large,  the  approximation  of  derivatives  by  secants  may  be 
very  poor  unless  the  particular  region  is  dlsected  into  smaller  regions  which  have 
reasonable  secant  approximations,  a practice  commonly  known  as  mesh  refinement.  In 
fluid  mechanics,  the  boundary  layer  of  a viscous  flow  around  or  through  an  object  is 
such  a region. 

Obviously,  the  necessary  resoltuion  could  be  accomplished  by  merely  increasing 
the  number  of  points  in  a uniform  distribution;  however,  this  would  require  excessive 
computer  time  and  storage.  Another  alternative,  known  as  the  interface  method,  is 
to  use  a refined  mesh  only  in  the  given  region  and  then  join  it  with  the  global 
mesh.  An  improved  technique  is  to  use  coordinate  distribution  functions  which 
smoothly  distribute  mesh  points  so  that  in  some  sense  they  are  spaced  in  roughly 
an  inverse  proportion  to  the  size  of  the  gradients.  Thus,  regions  of  high  gradients 
have  proportionately  more  points  than  regions  with  smaller  gradients.  Unlike  the 
interface  method,  the  transition  between  different  mesh  lengths  is  made  continuously, 
and  as  gradually  as  possible.  Distributions  are  often  used  when  the  distributional 
transformation  is  applied  to  an  independent  variable  of  an  existing  transformation. 

The  result  is  a new  transformation  obtained  by  composition.  With  this  approach,  the 
problem  of  mesh  point  distribution  is  replaced  by  the  problem  of  selecting  a suitable 
set  of  distribution  functions  within  a transformation  of  coordinates.  The  problem 
is  a nontrivial  one  since  the  distribution  functions  should  depend  upon  the  nature 
of  the  solution  being  computed  but  cure  determined  in  advance  of  the  computation. 

Thus,  some  prior  knowledge  of  the  solution  is  required . In  flows  with  large  boundary 
layer  separation  or  with  adjacent  dissimilar  components,  the  critical  region  to  be 
resolved  is  somewhere  in  the  middle  of  the  flow.  But  the  location  of  such,  regions 
is  often  unknown  at  the  outset  of  the  problem.  One  method  to  overcome  this  difficulty 
in  marching  procedures  is  to  create  the  distribution  function  at  the  next  level  based 
upon  a knowledge  of  the  solution  at  the  present  level.  Care  must  be  taken,  however, 
to  create  a distribution  function  that  is  sufficiently  smooth  in  the  marching 
direction. 

m many  problems  of  practical  interest,  however,  the  regions  that  need  resolu- 
tion are  known  in  advance.  Typical  examples  are  attached  boundary  layers  and  boundary 
layers  that  may  have  small  separations  or  separation  bubbles. 

Within  the  framework  of  tube-like  coordinate  Bystems  boundary  layer  resolution 
(Ref.  28)  on  the  inner  surface  is  accomplished  by  setting 
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where  d » d(t)  is  the  estimated  boundary  layer  growth,  or  « cr(t)  is  the  desired 
proportion  of  mesh  points  in  the  bound ary  layer,  and  D is  the  hyperbolic  damping 
factor.  The  boundary  layer  growth  d gives  the  fraction  of  the  flow  region  occupied 
by  the  boundary  layer,  a is  usually  taken  as  a constant,  and  D can  be  given  a value 
of  about  2.  When  r is  small,  the  radial  distribution  of  equation  (45)  reduces 
essentially  to  the  line 

R-  £ r (l6) 


which  would  have  been  chosen  had  we  used  the  interface  method.  As  r approaches 
unity  the  distribution  Bq.  (15)  smoothly  approaches  unity  as  illustrated  in  Fig.  6. 
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TREATMENT  OF  THE  BOW  SHOCK 


Although  the  shape  of  the  body  is  specified  at  the  start  of  the  problem,  the 
shape  of  the  bow  shock  cannot  be  specified  until  effects  of  the  body  on  the  flow 
within  the  bow  shock  are  evaluated.  The  shock  shape  is  therefore  evaluated  at  each 
new  marching  station  based  on  information  already  computed.  The  shock  surface  is 
adopted  as  the  outer  coordinate  surface  and  is  used  to  determine  the  necessary 
metric  information  for  the  tube-like  coordinates.  The  governing  equations  are  then 
solved  in  the  annular  region  between  the  ogival  body  and  the  bow  shock  by  marching 
from  one  transverse  plane  to  the  next,  proceeding  in  the  nominally  streamwise 
direction. 


The  bow  shock  is  computed  as  a discontinuity  satisfying  the  classical  Rankine- 
Hugoniot  relations.  The  intersection  of  this  shock  and  a transverse  computational 
plane  is  a loop  represented  by  discrete  grid  points.  Provided  that  a given  grid 
point  on  this  loop  is  outside  the  "zone  of  influence"  of  the  neighboring  points  on 
the  loop,  the  shock  solution  at  the  given  point  is  independent  of  the  solution  at 
adjacent  points.  This  "zone  of  influence"  assumption  is  valid  over  a wide  range 
of  flow  conditions  and  consequently  is  not  a limiting  assumption.  Thus  the  shock 
radius  yn+1(0 ) at  each  point  in  the  n+1  transverse  plane  can  be  evaluated  inde- 
pendently by  a pointwise  iteration  procedure. 


The  iteration  at  each  circumferential  location  in  the  n+1  plane  proceeds  by 
first  locally  extending  the  shock  surface  from  the  most  recently  evaluated  trans- 
verse computational  plane,  n,  to  a point  in  the  n+1  plane.  The  extension  of  the 
shock  surface  includes  the  point  being  evaluated  in  the  n+1  plane,  but  does  not 
extend  circumferentially  to  the  neighboring  points.  This  extension  is  a first 
guess  for  the  shock  location  at  a circumferential  point  and  hence  for  a point  on 
the  outer  tube-like  coordinate  surface  given  by  Yn+1(0)  in  Eq. (lUX 

Given  a guess  at  the  shock  location,  the  axial  mass  flux  inside  the  shock  can 
be  computed  by  two  methods.  First,  an  application  of  the  Rankine-Hugoniot  condi- 
tions produces  a value  of  the  axial  mass  flux  based  only  on  the  shock  shape  and  the 
flow  properties  outside  the  shock.  Second,  an  application  of  a compatibility  condi- 
tion produces  a second  value  of  this  flux  that  depends  only  on  the  shock  shape  and 
the  flow  properties  inside  the  shock.  The  shock  location  is  then  adjusted  iteratively 
until  the  axial  mass  flux  inside  the  shock  computed  by  the  two  methods  is  the  same. 
This  iteration  for  the  shock  location  is  repeated  at  each  of  the  circumferential 
grid  points  to  produce  a ring  of  discrete  points  at  the  n+1  station  which  collectively 
determine  the  shock  surface.  These  discrete  points  must  be  represented  by  a con- 
tinuous smooth  curve  to  provide  the  information  required  to  construct  the  coordinate 
system.  For  this  purpose  a least  squares-spline  curve  fitting  routine  is  employed. 
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RESULTS 


A number  of  test  cases  were  run  to  test,  evaluate,  and  extend  the  PEPSIG 
code  on  sequence  of  flow  problems.  The  sequence  of  problems  was  chosen  to 
reflect  the  various  types  of  flow  situations  that  are  likely  to  occur  in  the  com- 
plex flow  field  over  an  ogival  body  of  arbitrary  gjometry  when  it  is  at  an 
tingle  of  attack.  Since  the  flow  has  subsonic  regions,  supersonic  regions,  and  shock 
waves  it  is  imperative  that  test  cases  should  be  successfully  executed  when  these 
flow  types  are  isolated  and  when  they  occur  simultaneously.  Since  purely  subsonic 
flow  has  been  successfully  done  for  boundary  layers  and  for  internal  duct  flow 
under  a separate  project  for  NASA  Lewis,  it  need  not  be  displayed  here.  For 
purely  supersonic  flow  with  shocks,  several  test  cases  were  run.  First,  it  was 
desired  to  correctly  model  a shock  wave  and  allow  it  to  properly  be  passed  out  of 
the  computational  region.  The  ability  to  pass  shock  waves  out  of  a region  required 
the  application  of  exit  derivative  conditions  aligned  with  the  shock  angle.  In 
Fig.  7 a spurious  shock  reflection  is  observed  which  is  corrected  by  the  develop- 
ment of  exit  derivative  conditions  as  displayed  in  Fig.  8.  In  each  case,  the 
specified  shock  angle  was  correctly  preserved  as  the  solution  was  marched  upstream. 
However,  the  solution  in  each  case  possessed  wiggles.  To  remove  these  wiggles  a 
fourth  order  damping  term  was  added.  In  the  present  implicit  procedure  a modified 
form  of  the  fourth-order  pressure  damping  term  suggested  by  MacCormack  and  Baldwin 
(Ref.  31*)  and  modified  in  Ref.  35,  has  been  employed,  viz.,  a term  of  the  form 


D 


- „ Jit 

Pk  ^kdX,2 


(17) 


has  been  added  to  each  of  the  governing  equations  being  solved  for  each  coordinate 
direction  xk.  The  diffusion  coefficient  in  Eq.  (17 ) was  taken  as 


0(AXk)3f 


kk  c , 

1 

K2  / 

(18) 


where  0 < B - 1/2,  c is  the  local  sound  speed,  and  the  average  pressure  p given  by 

(«!_,♦  <»> 

where  i is  the  grid  point  index  in  the  x. -direction.  In  the  damping  term  the  factor 
f is  set  to  1.0  for  the  continuity  equation  and  to  p for  all  other  conservation 
equations.  For  the  free  shock,  cross  sectioned  results  without  damping  in  Figs.  9» 
and  10  are  improved  with  damping  as  displayed  in  Figs.  11  and  12.  Here,  a sequence 
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of  cross  sectional  cuts  are  given  to  show  that  the  remaining  post  shock  wiggles 
are  dissipated  into  smooth  curves  as  the  solution  is  marched  in  the  upstream 
direction.  With  the  ability  to  correctly  model  purely  subsonic  flow  and  purely 
supersonic  flow  with  shocks,  it  remains  to  demonstrate  test  cases  with  mixed 
subsonic  and  supersonic  regions.  The  initial  extension  in  this  direction  was 
performed  on  the  model  problem  of  a supersonic  flow  over  a flat  plate  with  no-slip 
boundary  conditions.  The  no-slip  boundary  condition  forces  a subsonic  region  in  the 
wall  vicinity.  In  particular,  a Mach  2 free  stream  flow  developed  a boundary  layer 
with  a subsonic  region.  Profiles  of  velocity,  density,  pressure,  and  Mach  number 
are  given  in  Figs.  13  to  17.  Figure  13  gives  calculated  values  of  the  streamwise 
velocity.  It  should  be  noted  that  although  a subsonic  region  clearly  exists  in  the 
vicinity  of  the  wall  (see  Fig.  17 ),  the  calculation  produced  qualitatively  reasonable 
profiles.  The  calculated  transverse  velocity,  density,  pressure  distributions  and 
Mach  number  are  presented  in  Figs.  lU-17,  respectively.  This  calculation  is  an 
indicator  of  a basic  capability  of  the  code  to  model  a supersonic  flow  with  a no-slip 
boundary  condition,  Similarly,  supersonic  flow  with  an  imbedded  shock  wave  over  a 
wedge  was  considered.  Cross  sectional  profiles  are  displayed  in  Figs.  18  and  19. 

In  each  case  above,  the  cross  sectional  profiles  appear  to  correctly  model  the  physics. 
However,  in  the  streamwise  direction,  there  was  a gradual  development  of  an  adverse 
pressure  gradient.  This  was  caused  by  supersonic  pressure  waves  impinging  upon  the 
subsonic  portion  (Refs.  25-26);  the  result  is  the  development  of  departure  solutions 
(Refs.  l6-l8,  25-26).  In  Issa  and  Lockwood  (Ref.  26)  the  suggestion  is  to  apply  a 
special  condition  at  the  sonic  line  while  the  other  investigators  (Refs.  16-18)  have 
tried  various  specifications  of  the  subsonic  axial  pressure  gradient.  With  due 
consideration  to  the  results  of  the  other  investigators,  a more  satisfactory 
resolution  to  the  problem  of  departure  solutions  is  still  being  sought. 
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